<!DOCTYPE html>

<html xmlns="http://www.w3.org/1999/xhtml">

<head>

<meta charset="utf-8" />
<meta http-equiv="Content-Type" content="text/html; charset=utf-8" />
<meta name="generator" content="pandoc" />

<meta name="viewport" content="width=device-width, initial-scale=1">

<meta name="author" content="Arezoo Rafieeinasab &amp; Aubrey Dugger" />

<meta name="date" content="2017-05-02" />

<title>Spatial tools</title>



<style type="text/css">code{white-space: pre;}</style>
<style type="text/css">
div.sourceCode { overflow-x: auto; }
table.sourceCode, tr.sourceCode, td.lineNumbers, td.sourceCode {
  margin: 0; padding: 0; vertical-align: baseline; border: none; }
table.sourceCode { width: 100%; line-height: 100%; }
td.lineNumbers { text-align: right; padding-right: 4px; padding-left: 4px; color: #aaaaaa; border-right: 1px solid #aaaaaa; }
td.sourceCode { padding-left: 5px; }
code > span.kw { color: #007020; font-weight: bold; } /* Keyword */
code > span.dt { color: #902000; } /* DataType */
code > span.dv { color: #40a070; } /* DecVal */
code > span.bn { color: #40a070; } /* BaseN */
code > span.fl { color: #40a070; } /* Float */
code > span.ch { color: #4070a0; } /* Char */
code > span.st { color: #4070a0; } /* String */
code > span.co { color: #60a0b0; font-style: italic; } /* Comment */
code > span.ot { color: #007020; } /* Other */
code > span.al { color: #ff0000; font-weight: bold; } /* Alert */
code > span.fu { color: #06287e; } /* Function */
code > span.er { color: #ff0000; font-weight: bold; } /* Error */
code > span.wa { color: #60a0b0; font-weight: bold; font-style: italic; } /* Warning */
code > span.cn { color: #880000; } /* Constant */
code > span.sc { color: #4070a0; } /* SpecialChar */
code > span.vs { color: #4070a0; } /* VerbatimString */
code > span.ss { color: #bb6688; } /* SpecialString */
code > span.im { } /* Import */
code > span.va { color: #19177c; } /* Variable */
code > span.cf { color: #007020; font-weight: bold; } /* ControlFlow */
code > span.op { color: #666666; } /* Operator */
code > span.bu { } /* BuiltIn */
code > span.ex { } /* Extension */
code > span.pp { color: #bc7a00; } /* Preprocessor */
code > span.at { color: #7d9029; } /* Attribute */
code > span.do { color: #ba2121; font-style: italic; } /* Documentation */
code > span.an { color: #60a0b0; font-weight: bold; font-style: italic; } /* Annotation */
code > span.cv { color: #60a0b0; font-weight: bold; font-style: italic; } /* CommentVar */
code > span.in { color: #60a0b0; font-weight: bold; font-style: italic; } /* Information */
</style>



<link href="data:text/css;charset=utf-8,body%20%7B%0Abackground%2Dcolor%3A%20%23fff%3B%0Amargin%3A%201em%20auto%3B%0Amax%2Dwidth%3A%20700px%3B%0Aoverflow%3A%20visible%3B%0Apadding%2Dleft%3A%202em%3B%0Apadding%2Dright%3A%202em%3B%0Afont%2Dfamily%3A%20%22Open%20Sans%22%2C%20%22Helvetica%20Neue%22%2C%20Helvetica%2C%20Arial%2C%20sans%2Dserif%3B%0Afont%2Dsize%3A%2014px%3B%0Aline%2Dheight%3A%201%2E35%3B%0A%7D%0A%23header%20%7B%0Atext%2Dalign%3A%20center%3B%0A%7D%0A%23TOC%20%7B%0Aclear%3A%20both%3B%0Amargin%3A%200%200%2010px%2010px%3B%0Apadding%3A%204px%3B%0Awidth%3A%20400px%3B%0Aborder%3A%201px%20solid%20%23CCCCCC%3B%0Aborder%2Dradius%3A%205px%3B%0Abackground%2Dcolor%3A%20%23f6f6f6%3B%0Afont%2Dsize%3A%2013px%3B%0Aline%2Dheight%3A%201%2E3%3B%0A%7D%0A%23TOC%20%2Etoctitle%20%7B%0Afont%2Dweight%3A%20bold%3B%0Afont%2Dsize%3A%2015px%3B%0Amargin%2Dleft%3A%205px%3B%0A%7D%0A%23TOC%20ul%20%7B%0Apadding%2Dleft%3A%2040px%3B%0Amargin%2Dleft%3A%20%2D1%2E5em%3B%0Amargin%2Dtop%3A%205px%3B%0Amargin%2Dbottom%3A%205px%3B%0A%7D%0A%23TOC%20ul%20ul%20%7B%0Amargin%2Dleft%3A%20%2D2em%3B%0A%7D%0A%23TOC%20li%20%7B%0Aline%2Dheight%3A%2016px%3B%0A%7D%0Atable%20%7B%0Amargin%3A%201em%20auto%3B%0Aborder%2Dwidth%3A%201px%3B%0Aborder%2Dcolor%3A%20%23DDDDDD%3B%0Aborder%2Dstyle%3A%20outset%3B%0Aborder%2Dcollapse%3A%20collapse%3B%0A%7D%0Atable%20th%20%7B%0Aborder%2Dwidth%3A%202px%3B%0Apadding%3A%205px%3B%0Aborder%2Dstyle%3A%20inset%3B%0A%7D%0Atable%20td%20%7B%0Aborder%2Dwidth%3A%201px%3B%0Aborder%2Dstyle%3A%20inset%3B%0Aline%2Dheight%3A%2018px%3B%0Apadding%3A%205px%205px%3B%0A%7D%0Atable%2C%20table%20th%2C%20table%20td%20%7B%0Aborder%2Dleft%2Dstyle%3A%20none%3B%0Aborder%2Dright%2Dstyle%3A%20none%3B%0A%7D%0Atable%20thead%2C%20table%20tr%2Eeven%20%7B%0Abackground%2Dcolor%3A%20%23f7f7f7%3B%0A%7D%0Ap%20%7B%0Amargin%3A%200%2E5em%200%3B%0A%7D%0Ablockquote%20%7B%0Abackground%2Dcolor%3A%20%23f6f6f6%3B%0Apadding%3A%200%2E25em%200%2E75em%3B%0A%7D%0Ahr%20%7B%0Aborder%2Dstyle%3A%20solid%3B%0Aborder%3A%20none%3B%0Aborder%2Dtop%3A%201px%20solid%20%23777%3B%0Amargin%3A%2028px%200%3B%0A%7D%0Adl%20%7B%0Amargin%2Dleft%3A%200%3B%0A%7D%0Adl%20dd%20%7B%0Amargin%2Dbottom%3A%2013px%3B%0Amargin%2Dleft%3A%2013px%3B%0A%7D%0Adl%20dt%20%7B%0Afont%2Dweight%3A%20bold%3B%0A%7D%0Aul%20%7B%0Amargin%2Dtop%3A%200%3B%0A%7D%0Aul%20li%20%7B%0Alist%2Dstyle%3A%20circle%20outside%3B%0A%7D%0Aul%20ul%20%7B%0Amargin%2Dbottom%3A%200%3B%0A%7D%0Apre%2C%20code%20%7B%0Abackground%2Dcolor%3A%20%23f7f7f7%3B%0Aborder%2Dradius%3A%203px%3B%0Acolor%3A%20%23333%3B%0Awhite%2Dspace%3A%20pre%2Dwrap%3B%20%0A%7D%0Apre%20%7B%0Aborder%2Dradius%3A%203px%3B%0Amargin%3A%205px%200px%2010px%200px%3B%0Apadding%3A%2010px%3B%0A%7D%0Apre%3Anot%28%5Bclass%5D%29%20%7B%0Abackground%2Dcolor%3A%20%23f7f7f7%3B%0A%7D%0Acode%20%7B%0Afont%2Dfamily%3A%20Consolas%2C%20Monaco%2C%20%27Courier%20New%27%2C%20monospace%3B%0Afont%2Dsize%3A%2085%25%3B%0A%7D%0Ap%20%3E%20code%2C%20li%20%3E%20code%20%7B%0Apadding%3A%202px%200px%3B%0A%7D%0Adiv%2Efigure%20%7B%0Atext%2Dalign%3A%20center%3B%0A%7D%0Aimg%20%7B%0Abackground%2Dcolor%3A%20%23FFFFFF%3B%0Apadding%3A%202px%3B%0Aborder%3A%201px%20solid%20%23DDDDDD%3B%0Aborder%2Dradius%3A%203px%3B%0Aborder%3A%201px%20solid%20%23CCCCCC%3B%0Amargin%3A%200%205px%3B%0A%7D%0Ah1%20%7B%0Amargin%2Dtop%3A%200%3B%0Afont%2Dsize%3A%2035px%3B%0Aline%2Dheight%3A%2040px%3B%0A%7D%0Ah2%20%7B%0Aborder%2Dbottom%3A%204px%20solid%20%23f7f7f7%3B%0Apadding%2Dtop%3A%2010px%3B%0Apadding%2Dbottom%3A%202px%3B%0Afont%2Dsize%3A%20145%25%3B%0A%7D%0Ah3%20%7B%0Aborder%2Dbottom%3A%202px%20solid%20%23f7f7f7%3B%0Apadding%2Dtop%3A%2010px%3B%0Afont%2Dsize%3A%20120%25%3B%0A%7D%0Ah4%20%7B%0Aborder%2Dbottom%3A%201px%20solid%20%23f7f7f7%3B%0Amargin%2Dleft%3A%208px%3B%0Afont%2Dsize%3A%20105%25%3B%0A%7D%0Ah5%2C%20h6%20%7B%0Aborder%2Dbottom%3A%201px%20solid%20%23ccc%3B%0Afont%2Dsize%3A%20105%25%3B%0A%7D%0Aa%20%7B%0Acolor%3A%20%230033dd%3B%0Atext%2Ddecoration%3A%20none%3B%0A%7D%0Aa%3Ahover%20%7B%0Acolor%3A%20%236666ff%3B%20%7D%0Aa%3Avisited%20%7B%0Acolor%3A%20%23800080%3B%20%7D%0Aa%3Avisited%3Ahover%20%7B%0Acolor%3A%20%23BB00BB%3B%20%7D%0Aa%5Bhref%5E%3D%22http%3A%22%5D%20%7B%0Atext%2Ddecoration%3A%20underline%3B%20%7D%0Aa%5Bhref%5E%3D%22https%3A%22%5D%20%7B%0Atext%2Ddecoration%3A%20underline%3B%20%7D%0A%0Acode%20%3E%20span%2Ekw%20%7B%20color%3A%20%23555%3B%20font%2Dweight%3A%20bold%3B%20%7D%20%0Acode%20%3E%20span%2Edt%20%7B%20color%3A%20%23902000%3B%20%7D%20%0Acode%20%3E%20span%2Edv%20%7B%20color%3A%20%2340a070%3B%20%7D%20%0Acode%20%3E%20span%2Ebn%20%7B%20color%3A%20%23d14%3B%20%7D%20%0Acode%20%3E%20span%2Efl%20%7B%20color%3A%20%23d14%3B%20%7D%20%0Acode%20%3E%20span%2Ech%20%7B%20color%3A%20%23d14%3B%20%7D%20%0Acode%20%3E%20span%2Est%20%7B%20color%3A%20%23d14%3B%20%7D%20%0Acode%20%3E%20span%2Eco%20%7B%20color%3A%20%23888888%3B%20font%2Dstyle%3A%20italic%3B%20%7D%20%0Acode%20%3E%20span%2Eot%20%7B%20color%3A%20%23007020%3B%20%7D%20%0Acode%20%3E%20span%2Eal%20%7B%20color%3A%20%23ff0000%3B%20font%2Dweight%3A%20bold%3B%20%7D%20%0Acode%20%3E%20span%2Efu%20%7B%20color%3A%20%23900%3B%20font%2Dweight%3A%20bold%3B%20%7D%20%20code%20%3E%20span%2Eer%20%7B%20color%3A%20%23a61717%3B%20background%2Dcolor%3A%20%23e3d2d2%3B%20%7D%20%0A" rel="stylesheet" type="text/css" />

</head>

<body>




<h1 class="title toc-ignore">Spatial tools</h1>
<h4 class="author"><em>Arezoo Rafieeinasab &amp; Aubrey Dugger</em></h4>
<h4 class="date"><em>2017-05-02</em></h4>



<p>For model analysis and evaluation, we often need to create spatial maps, aggregate over spatial units, or produce georeferenced rasters and shapefiles. We have adapted existing functionality from spatial libraries such as SP, RGDAL, RGEOS, and Raster into rwrfhydro. In this vignette, we describe some of these spatial functions and give examples of their application to model evaluation.</p>
<div id="list-of-the-available-functions" class="section level2">
<h2>List of the available functions</h2>
<ul>
<li>GetProj</li>
<li>GetGeogridSpatialInfo</li>
<li>ExportGeogrid</li>
<li>GetGeogridIndex</li>
<li>GetTimeZone</li>
<li>GetRfc</li>
<li>GetPoly</li>
<li>PolygonToRaster</li>
</ul>
</div>
<div id="general-info" class="section level2">
<h2>General Info</h2>
<p>Load the rwrfhydro package.</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="kw">library</span>(rwrfhydro)
<span class="kw">library</span>(rgdal) ## on some linux machines this appears needed
<span class="kw">options</span>(<span class="dt">warn=</span><span class="dv">1</span>)</code></pre></div>
<p>Set a data path to the Fourmile Creek test case.</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r">fcPath &lt;-<span class="st"> '~/wrfHydroTestCases/Fourmile_Creek_testcase_v2.0'</span></code></pre></div>
<p>The geogrid file is the main coarse-grid (LSM) parameter file and contains base geographic information on the model domain such as the geographic coordinate system and latitude/longitude coordinates. We use this file frequently. Set a path to geogrid file.</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r">geoFile &lt;-<span class="st"> </span><span class="kw">paste0</span>(fcPath,<span class="st">'/DOMAIN/geo_em_d01.Fourmile1km.nlcd11.nc'</span>)</code></pre></div>
</div>
<div id="getproj" class="section level2">
<h2>GetProj</h2>
<p>To be able to use spatial tools in R, we need to know the projection information for the domain. All the coarse-resolution (LSM) model input and output files are based on the geogrid domain. <code>GetProj</code> will pull projection information from the geogrid file. You will see the specs for the Lambert Conformal Conic projection on the WRF standard sphere.</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r">proj4 &lt;-<span class="st"> </span><span class="kw">GetProj</span>(geoFile)
proj4</code></pre></div>
<pre><code>## [1] &quot;+proj=lcc +lat_1=30 +lat_2=50 +lat_0=40.0403709411621 +lon_0=-105 +x_0=0 +y_0=0 +a=6370000 +b=6370000 +units=m +no_defs&quot;</code></pre>
</div>
<div id="getgeogridspatialinfo" class="section level2">
<h2>GetGeogridSpatialInfo</h2>
<p><code>GetGeogridSpatialInfo</code> will pull geospatial information about the coarse-resolution (LSM) model domain from the geogrid file.</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r">geoInfo &lt;-<span class="st"> </span><span class="kw">GetGeogridSpatialInfo</span>(geoFile)
geoInfo</code></pre></div>
<pre><code>##     DX   DY GRID_TYPE     LAT1      LON1  REF_LAT REF_LON SPLAT1 SPLAT2
## 1 1000 1000         C 40.01251 -105.5616 40.04037    -105     30     50
##   POLE_LAT POLE_LON MAP_PROJ NCOL NROW
## 1       90        0        1   21    7</code></pre>
</div>
<div id="exportgeogrid" class="section level2">
<h2>ExportGeogrid</h2>
<p>If you need to create a georeferenced TIF file from any variable in an LSM-related netcdf file (input or output), then you can use the <code>ExportGeogrid</code> function. It takes a NetCDF file and converts the specified variable into a georeferenced TIF file for use in standard GIS tools. You can use <code>ExportGeogrid</code> directly on a file that contains lat/lon coordinates or you can use it on a file that does not contain lat/lon coords by providing a separate coordinate file.</p>
<p>Let’s export a variable from the geogrid file. You can get a list of all available variables in the <code>geoFile</code> using the <code>ncdump</code> function in rwrfhydro.</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="kw">head</span>(<span class="kw">ncdump</span>(geoFile))</code></pre></div>
<p>Now we will create a georeferenced TIF file from the elevation field. The geogrid contains lat/long coordinates, so you only need to provide the address to the file (<code>geoFile</code>), the name of the variable (<code>HGT_M</code>), and the name of the output file (<code>geogrid_hgt.tif</code>).</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="kw">ExportGeogrid</span>(geoFile,<span class="st">&quot;HGT_M&quot;</span>, <span class="st">&quot;geogrid_hgt.tif&quot;</span>)</code></pre></div>
<p>You can now use the geotiff in any standard GIS platform. Here we will just read it into memory as a raster and display it.</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="co"># Read the newly created tiff file</span>
<span class="kw">library</span>(raster)
r &lt;-<span class="st"> </span><span class="kw">raster</span>(<span class="st">&quot;geogrid_hgt.tif&quot;</span>)

<span class="co"># Plot the imported raster from tiff file</span>
<span class="kw">plot</span>(r, <span class="dt">main =</span> <span class="st">&quot;HGT_M&quot;</span>, <span class="dt">col=</span><span class="kw">terrain.colors</span>(<span class="dv">100</span>))

<span class="co"># Check the raser information and notice that geographic coordinate information has been added.</span>
r</code></pre></div>
<pre><code>## class       : RasterLayer 
## dimensions  : 7, 21, 147  (nrow, ncol, ncell)
## resolution  : 1000, 1000  (x, y)
## extent      : -47596.92, -26596.92, -3400.022, 3599.978  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=lcc +lat_1=30 +lat_2=50 +lat_0=40.0403709411621 +lon_0=-105 +x_0=0 +y_0=0 +a=6370000 +b=6370000 +units=m +no_defs 
## data source : /Users/james/R/jlm_lib/rwrfhydro/vignettes/geogrid_hgt.tif 
## names       : geogrid_hgt</code></pre>
<p><img src="" width="600" height="600" /></p>
<p>Many of the output files (such as LDASOUT, RESTARTS) do not contain lat/lon coordinates but match the spatial coordinate system of the geogrid input file. In that case, you can provide a supplemental <code>inCoordFile</code> which contains the lat/lon information.</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r">file &lt;-<span class="st"> </span><span class="kw">paste0</span>(fcPath,<span class="st">&quot;/run.FluxEval/RESTART.2013060100_DOMAIN1&quot;</span>)
<span class="co"># ncdump(file) # check if the SOIL_T exist in the file</span>

<span class="co"># Export the 3rd layer of the 4-layer soil temperature variable</span>
<span class="kw">ExportGeogrid</span>(file,
             <span class="dt">inVar=</span><span class="st">&quot;SOIL_T&quot;</span>,
             <span class="dt">outFile=</span><span class="st">&quot;20130315_soilm3.tif&quot;</span>,
             <span class="dt">inCoordFile=</span>geoFile,
             <span class="dt">inLyr=</span><span class="dv">3</span>)

<span class="co"># Read the newly created tiff file</span>
r &lt;-<span class="st"> </span><span class="kw">raster</span>(<span class="st">&quot;20130315_soilm3.tif&quot;</span>)

<span class="co"># Plot the imported raster from tiff file</span>
<span class="kw">plot</span>(r, <span class="dt">main =</span> <span class="st">&quot;Soil Temperature&quot;</span>, <span class="dt">col=</span><span class="kw">rev</span>(<span class="kw">heat.colors</span>(<span class="dv">100</span>))) <span class="co"># in raster</span></code></pre></div>
<p><img src="" width="600" height="600" /></p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="co"># Check the raster information and notice that geographic coordinate information has been added</span>
r</code></pre></div>
<pre><code>## class       : RasterLayer 
## dimensions  : 7, 21, 147  (nrow, ncol, ncell)
## resolution  : 1000, 1000  (x, y)
## extent      : -47596.92, -26596.92, -3400.022, 3599.978  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=lcc +lat_1=30 +lat_2=50 +lat_0=40.0403709411621 +lon_0=-105 +x_0=0 +y_0=0 +a=6370000 +b=6370000 +units=m +no_defs 
## data source : /Users/james/R/jlm_lib/rwrfhydro/vignettes/20130315_soilm3.tif 
## names       : X20130315_soilm3</code></pre>
</div>
<div id="getgeogridindex" class="section level2">
<h2>GetGeogridIndex</h2>
<p>To be able to use tools such as <code>GetMultiNcdf</code> to pull data from gridded output, we need to know the indices (i,j) of the area of interest within the domain. <code>GetGeogridIndex</code> calculates cell indices from lat/lon (or other) coordinates. It reads in a set of lat/lon (or other) coordinates and generates a corresponding set of geogrid index pairs. You can assign a projection to the points using the <code>proj4</code> argument, which will be used to transform the points to the <code>geoFile</code> coordinate system. Check <code>?GetGeogridIndex</code> or the Precipitation Evaluation vignette for full usage.</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r">sg &lt;-<span class="st"> </span><span class="kw">data.frame</span>(<span class="dt">lon =</span> <span class="kw">seq</span>(<span class="op">-</span><span class="fl">105.562</span>, <span class="op">-</span><span class="fl">105.323</span>, <span class="dt">length.out =</span> <span class="dv">10</span>), 
                 <span class="dt">lat =</span> <span class="kw">seq</span>(<span class="fl">40.0125</span>, <span class="fl">40.0682</span>, <span class="dt">length.out =</span> <span class="dv">10</span>))
<span class="kw">GetGeogridIndex</span>(sg, geoFile)</code></pre></div>
<pre><code>##    we sn
## 1   1  1
## 2   3  2
## 3   5  2
## 4   8  3
## 5  10  4
## 6  12  4
## 7  14  5
## 8  17  6
## 9  19  6
## 10 21  7</code></pre>
</div>
<div id="gettimezone" class="section level2">
<h2>GetTimeZone</h2>
<p>Many station observations are reported in local time and need to be converted to UTC time to be comparable with WRF-Hydro inputs and outputs. <code>GetTimeZone</code> returns the time zone for any lat/lon coordinates. It simply takes a dataframe containing at least two fields of <code>latitude</code> and <code>longitude</code> and overlays the <code>points</code> with a timezone shapefile (can be downloded from <a href="http://efele.net/maps/tz/world/" class="uri">http://efele.net/maps/tz/world/</a>). The shapefile is provided in rwrfhydro data and is called <code>timeZone</code>.</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="co"># timeZone has been provided by rwrfhydro as a SpatialPolygonDataFrame</span>
<span class="kw">class</span>(timeZone)</code></pre></div>
<pre><code>## [1] &quot;SpatialPolygonsDataFrame&quot;
## attr(,&quot;package&quot;)
## [1] &quot;sp&quot;</code></pre>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="co"># Shows the available timezone (TZID column in timeZone@data)</span>
<span class="kw">head</span>(timeZone<span class="op">@</span>data)</code></pre></div>
<pre><code>##                  TZID
## 0    America/Dominica
## 1  America/St_Vincent
## 2 Australia/Lord_Howe
## 3        Asia/Kashgar
## 4      Pacific/Wallis
## 5      Asia/Jerusalem</code></pre>
<p><code>GetTimeZone</code> has three arguments.</p>
<ul>
<li><code>points</code>: A dataframe of the points. The dataframe should contain at least two fields called <code>latitude</code> and <code>longitude</code>.</li>
<li><code>proj4</code>: Projection of the <code>points</code> to be used in transforming the <code>points</code> projection to the <code>timeZone</code> projection. Default is <code>+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0</code> which is the same as the <code>timezone</code> projection.</li>
<li><code>parallel</code>: If the number of points is high you can parallelize the process.</li>
</ul>
<p><code>GetTimeZone</code> will return the <code>points</code> dataframe with an added column called <code>timeZone</code>. It will return NA if a point is not in any polygon. Now let’s generate some points and find their time zone information.</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="co"># Provide a dataframe of 10 points having longitude and latitude as column name.</span>
sg &lt;-<span class="st"> </span><span class="kw">data.frame</span>(<span class="dt">longitude =</span> <span class="kw">seq</span>(<span class="op">-</span><span class="dv">110</span>, <span class="op">-</span><span class="dv">80</span>, <span class="dt">length.out =</span> <span class="dv">10</span>),
                 <span class="dt">latitude =</span> <span class="kw">seq</span>(<span class="dv">30</span>, <span class="dv">50</span>, <span class="dt">length.out =</span> <span class="dv">10</span>))

<span class="co"># Find the time zone for each point</span>
sg &lt;-<span class="st"> </span><span class="kw">GetTimeZone</span>(sg)
sg</code></pre></div>
<pre><code>##              timeZone  longitude latitude
## 1  America/Hermosillo -110.00000 30.00000
## 2      America/Denver -106.66667 32.22222
## 3      America/Denver -103.33333 34.44444
## 4     America/Chicago -100.00000 36.66667
## 5     America/Chicago  -96.66667 38.88889
## 6     America/Chicago  -93.33333 41.11111
## 7     America/Chicago  -90.00000 43.33333
## 8                &lt;NA&gt;  -86.66667 45.55556
## 9     America/Toronto  -83.33333 47.77778
## 10    America/Toronto  -80.00000 50.00000</code></pre>
</div>
<div id="getrfc" class="section level2">
<h2>GetRfc</h2>
<p>The US has 13 River Forecast Centers (RFCs) which issue daily river forecasts using hydrologic models based on rainfall, soil characteristics, precipitation forecasts, and several other variables. The RFC boundary shapefile is provided in rwrfhydro data and is called <code>rfc</code>.</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="kw">class</span>(rfc)</code></pre></div>
<pre><code>## [1] &quot;SpatialPolygonsDataFrame&quot;
## attr(,&quot;package&quot;)
## [1] &quot;sp&quot;</code></pre>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="co"># Shows the available rfc, name of the column is BASIN_ID</span>
<span class="kw">head</span>(rfc<span class="op">@</span>data)</code></pre></div>
<pre><code>##   SITE_ID STATE           RFC_NAME                  RFC_CITY BASIN_ID
## 0     ACR    AK             Alaska                 Anchorage    AKRFC
## 1     KRF    MO     Missouri Basin Kansas City/Pleasant Hill    MBRFC
## 2     STR    UT     Colorado Basin            Salt Lake City    CBRFC
## 3     TUA    OK Arkansas-Red Basin                     Tulsa    ABRFC
## 4     RSA    CA  California-Nevada                Sacramento    CNRFC
## 5     ORN    LA  Lower Mississippi   New Orleans/Baton Rouge    LMRFC</code></pre>
<p><code>GetRfc</code> return the RFC name for any point having <code>latitude</code> and <code>longitude</code>. It takes a dataframe containing at least two fields of <code>latitude</code> and <code>longitude</code>, overlays the points with the <code>rfc</code> SpatialPolygonDataFrame, and return the RFC’s BASIN_ID. This function has three arguments:</p>
<ul>
<li><code>points</code>: A dataframe of the points. The dataframe should contain at least two fields called “latitude” and “longitude”.</li>
<li><code>proj4</code>: Projection of the <code>points</code> to be used in transforming the <code>points</code> projection to the <code>rfc</code> projection. Default is <code>+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0</code>.</li>
<li><code>parallel</code>: If the number of points is high you can parallelize the process.</li>
</ul>
<p><code>GetRfc</code> will return the points dataframe with an added column called <code>rfc</code>. It will return NA if the point is not in any polygon.</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="co"># Provide a dataframe of 10 points having longitude and latitude as column name.</span>
sg &lt;-<span class="st"> </span><span class="kw">data.frame</span>(<span class="dt">longitude =</span> <span class="kw">seq</span>(<span class="op">-</span><span class="dv">110</span>, <span class="op">-</span><span class="dv">80</span>, <span class="dt">length.out =</span> <span class="dv">10</span>), 
                 <span class="dt">latitude =</span> <span class="kw">seq</span>(<span class="dv">30</span>, <span class="dv">50</span>, <span class="dt">length.out =</span> <span class="dv">10</span>))

<span class="co"># Find the rfc for each point</span>
sg &lt;-<span class="st"> </span><span class="kw">GetRfc</span>(sg)
sg</code></pre></div>
<pre><code>##      rfc  longitude latitude
## 1   &lt;NA&gt; -110.00000 30.00000
## 2  WGRFC -106.66667 32.22222
## 3  WGRFC -103.33333 34.44444
## 4  ABRFC -100.00000 36.66667
## 5  MBRFC  -96.66667 38.88889
## 6  NCRFC  -93.33333 41.11111
## 7  NCRFC  -90.00000 43.33333
## 8   &lt;NA&gt;  -86.66667 45.55556
## 9   &lt;NA&gt;  -83.33333 47.77778
## 10  &lt;NA&gt;  -80.00000 50.00000</code></pre>
</div>
<div id="getpoly" class="section level2">
<h2>GetPoly</h2>
<p><code>GetPoly</code> is similar to <code>GetRfc</code>; it is a wrapper for the function <code>sp::over</code>. It takes a dataframe containing at least two fields of <code>latitude</code> and <code>longitude</code>, overlays the points with a <code>SpatialPolygonDataFrame</code>, and returns the requested attribute from the polygon. You could use any available <code>SpatialPolygon*</code> loaded into memory or provide the address to the location of a polygon shapefile and it will read the shapefile using the <code>rgdal::readOGR</code> function.</p>
<p>Let’s get the RFC information from <code>GetPoly</code> instead of <code>GetRfc</code>. Here we provide the name of the <code>SpatialPolygon*</code> and, using the argument <code>join</code>, request one of the polygon attributes. For example, here we request the <code>BASIN_ID</code>, <code>RFC_NAME</code> and <code>RFC_CITY</code> attributes.</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="co"># Provide a dataframe of 10 points having longitude and latitude</span>
sg &lt;-<span class="st"> </span><span class="kw">data.frame</span>(<span class="dt">longitude =</span> <span class="kw">seq</span>(<span class="op">-</span><span class="dv">110</span>, <span class="op">-</span><span class="dv">80</span>, <span class="dt">length.out =</span> <span class="dv">10</span>), 
                 <span class="dt">latitude =</span> <span class="kw">seq</span>(<span class="dv">30</span>, <span class="dv">50</span>, <span class="dt">length.out =</span> <span class="dv">10</span>))

<span class="co"># Find the ID of RFC for each point</span>
sg &lt;-<span class="st"> </span><span class="kw">GetPoly</span>(<span class="dt">points =</span> sg, <span class="dt">polygon =</span> rfc, <span class="dt">join =</span> <span class="st">&quot;BASIN_ID&quot;</span>)

<span class="co"># Find the full name of RFC for each point</span>
sg &lt;-<span class="st"> </span><span class="kw">GetPoly</span>(<span class="dt">points =</span> sg, <span class="dt">polygon =</span> rfc, <span class="dt">join =</span> <span class="st">&quot;RFC_NAME&quot;</span>)

<span class="co"># Find the location/city of RFC for each point</span>
sg &lt;-<span class="st"> </span><span class="kw">GetPoly</span>(<span class="dt">points =</span> sg, <span class="dt">polygon =</span> rfc, <span class="dt">join =</span> <span class="st">&quot;RFC_CITY&quot;</span>)
sg</code></pre></div>
<pre><code>##    BASIN_ID           RFC_NAME                  RFC_CITY  longitude
## 1      &lt;NA&gt;               &lt;NA&gt;                      &lt;NA&gt; -110.00000
## 2     WGRFC          West Gulf         Dallas/Fort Worth -106.66667
## 3     WGRFC          West Gulf         Dallas/Fort Worth -103.33333
## 4     ABRFC Arkansas-Red Basin                     Tulsa -100.00000
## 5     MBRFC     Missouri Basin Kansas City/Pleasant Hill  -96.66667
## 6     NCRFC      North Central               Minneapolis  -93.33333
## 7     NCRFC      North Central               Minneapolis  -90.00000
## 8      &lt;NA&gt;               &lt;NA&gt;                      &lt;NA&gt;  -86.66667
## 9      &lt;NA&gt;               &lt;NA&gt;                      &lt;NA&gt;  -83.33333
## 10     &lt;NA&gt;               &lt;NA&gt;                      &lt;NA&gt;  -80.00000
##    latitude
## 1  30.00000
## 2  32.22222
## 3  34.44444
## 4  36.66667
## 5  38.88889
## 6  41.11111
## 7  43.33333
## 8  45.55556
## 9  47.77778
## 10 50.00000</code></pre>
<p>Now let’s provide the address to a shapefile on the disk as well as the name of the shapefile and perform the same process. We have clipped the <code>HUC12</code> shapefile and provided it in the case study as a sample. The northeast portion of the clipped polygon partially covers the Fourmile Creek domain.</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="co"># Provide a dataframe of 10 points within the Fourmile Creek domain having longitude and latitude</span>
sg &lt;-<span class="st"> </span><span class="kw">data.frame</span>(<span class="dt">longitude =</span> <span class="kw">seq</span>(<span class="op">-</span><span class="fl">105.562</span>, <span class="op">-</span><span class="fl">105.323</span>, <span class="dt">length.out =</span> <span class="dv">10</span>), 
                 <span class="dt">latitude =</span> <span class="kw">seq</span>(<span class="fl">40.0125</span>, <span class="fl">40.0682</span>, <span class="dt">length.out =</span> <span class="dv">10</span>))


<span class="co"># We use `rgdal::readOG` in the GetPoly function and it does not interpret the character/symbol `~`. </span>
<span class="co"># Therefore, we need to use path.expand to get the full address to the case study location on your system. </span>
polygonAddress &lt;-<span class="st"> </span><span class="kw">paste0</span>(<span class="kw">path.expand</span>(fcPath), <span class="st">&quot;/polygons&quot;</span>)


<span class="co"># Find the HUC12 for each point</span>
sg &lt;-<span class="st"> </span><span class="kw">GetPoly</span>(<span class="dt">points =</span> sg,
              <span class="dt">polygonAddress =</span> polygonAddress,
              <span class="dt">polygonShapeFile =</span> <span class="st">&quot;clipped_huc12&quot;</span>,
              <span class="dt">join =</span> <span class="st">&quot;HUC12&quot;</span>)</code></pre></div>
<pre><code>## OGR data source with driver: ESRI Shapefile 
## Source: &quot;/Users/james/wrfHydroTestCases/Fourmile_Creek_testcase_v2.0/polygons&quot;, layer: &quot;clipped_huc12&quot;
## with 90 features
## It has 20 fields
## Integer64 fields read as strings:  OBJECTID GNIS_ID</code></pre>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r">sg</code></pre></div>
<pre><code>##           HUC12 longitude latitude
## 1  101900050401 -105.5620 40.01250
## 2  101900050401 -105.5354 40.01869
## 3  101900050403 -105.5089 40.02488
## 4  101900050403 -105.4823 40.03107
## 5  101900050403 -105.4558 40.03726
## 6  101900050403 -105.4292 40.04344
## 7  101900050403 -105.4027 40.04963
## 8  101900050403 -105.3761 40.05582
## 9          &lt;NA&gt; -105.3496 40.06201
## 10         &lt;NA&gt; -105.3230 40.06820</code></pre>
</div>
<div id="polytoraster" class="section level2">
<h2>PolyToRaster</h2>
<p>If you want to create an area mask in the coarse-resolution (LSM) model domain, you can use <code>PolyToRaster</code>. It first picks up the required geographic information (like <code>proj4</code>) from the geogrid file (<code>geoFile</code>) and then uses the <code>raster::rasterize</code> function to grab the mask or attibute values from the <code>SpatialPolygonDataFrame</code>. This function is basically wrapping the <code>raster::rasterize</code> function to serve our purpose. Below are a few different ways we can use this function.</p>
<p>Example 1 : Let’s get the RFC ID for each pixel within the Fourmile Creek domain. This is equivalent to rasterizing the <code>rfc</code> <code>SpatialPolygonDataFrame</code> based on the <code>BASIN_ID</code>.</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r">r &lt;-<span class="st"> </span><span class="kw">PolyToRaster</span>(<span class="dt">geoFile =</span> geoFile,
                  <span class="dt">useRfc =</span> <span class="ot">TRUE</span>,
                  <span class="dt">field =</span><span class="st">&quot;BASIN_ID&quot;</span>)</code></pre></div>
<pre><code>## Loading required namespace: rgeos</code></pre>
<p><img src="" width="600" height="600" /></p>
<p>To get the string (RFC ID) associated with the gridded ID value, you can check the attributes.</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r">r<span class="op">@</span>data<span class="op">@</span>attributes </code></pre></div>
<pre><code>## [[1]]
##   ID BASIN_ID
## 1  7    MBRFC</code></pre>
<p>As the results show, the full domain falls into one RFC (MBRFC).</p>
<p>Example 2 : Rasterize the HUC12 <code>SpatialPolygonDataFrame</code> based on the <code>HUC12</code> field. The clipped HUC12 shapefile is provided with the test case. You can read the shapefile and plot it as below.</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r">polyg &lt;-<span class="st"> </span>rgdal<span class="op">::</span><span class="kw">readOGR</span>(<span class="kw">paste0</span>(<span class="kw">path.expand</span>(fcPath), <span class="st">&quot;/polygons&quot;</span>), <span class="st">&quot;clipped_huc12&quot;</span>)
<span class="kw">plot</span>(polyg, <span class="dt">main =</span> <span class="st">&quot;Clipped HUC12&quot;</span>) ## in raster</code></pre></div>
<p><img src="" width="600" height="600" /></p>
<p>Our study domain partially covers a few basins in the northeast portion of this shapefile.</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r">polygonAddress &lt;-<span class="st"> </span><span class="kw">paste0</span>(<span class="kw">path.expand</span>(fcPath), <span class="st">&quot;/polygons&quot;</span>)
r &lt;-<span class="st"> </span><span class="kw">PolyToRaster</span>(<span class="dt">geoFile =</span> geoFile,
                  <span class="dt">polygonAddress =</span> polygonAddress,
                  <span class="dt">polygonShapeFile =</span> <span class="st">&quot;clipped_huc12&quot;</span>,
                  <span class="dt">field =</span><span class="st">&quot;HUC12&quot;</span>)</code></pre></div>
<p><img src="" width="600" height="600" /></p>
<p>To get the <code>HUC12</code> actual values:</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r">r<span class="op">@</span>data<span class="op">@</span>attributes</code></pre></div>
<pre><code>## [[1]]
##   ID        HUC12
## 1 33 101900050401
## 2 35 101900050403
## 3 36 101900050404
## 4 37 101900050406</code></pre>
<p>Example 3: You can create masks over the study domain using PolyToRaster. To create a unified mask over the study domain:</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r">r &lt;-<span class="st"> </span><span class="kw">PolyToRaster</span>(<span class="dt">geoFile =</span> geoFile,
                  <span class="dt">polygonAddress =</span> polygonAddress,
                  <span class="dt">polygonShapeFile =</span> <span class="st">&quot;clipped_huc12&quot;</span>,
                  <span class="dt">mask =</span><span class="ot">TRUE</span>)</code></pre></div>
<p><img src="" width="600" height="600" /></p>
<p>You can also create a separate mask for each subbasin (HUC12 in this case) with the fraction of each grid cell that is covered by each polygon. The fraction covered is estimated by dividing each cell into 100 subcells and determining presence/absence of the polygon in the center of each subcell.</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r">r &lt;-<span class="st"> </span><span class="kw">PolyToRaster</span>(<span class="dt">geoFile =</span> geoFile,
                  <span class="dt">polygonAddress =</span> polygonAddress,
                  <span class="dt">polygonShapeFile =</span> <span class="st">&quot;clipped_huc12&quot;</span>,
                  <span class="dt">field =</span> <span class="st">&quot;HUC12&quot;</span>,
                  <span class="dt">getCover =</span> <span class="ot">TRUE</span>)
<span class="kw">plot</span>(r) ## in raster</code></pre></div>
<p><img src="" width="600" height="600" /></p>
</div>



<!-- dynamically load mathjax for compatibility with self-contained -->
<script>
  (function () {
    var script = document.createElement("script");
    script.type = "text/javascript";
    script.src  = "https://mathjax.rstudio.com/latest/MathJax.js?config=TeX-AMS-MML_HTMLorMML";
    document.getElementsByTagName("head")[0].appendChild(script);
  })();
</script>

</body>
</html>
